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ABSTRACT 

We review some applications of fractional calculus developed by the author (partly in 
collaboration with others) to treat some basic problems in continuum and statistical 
mechanics. The problems in continuum mechanics concern mathematical modelling 
of viscoelastic bodies (§1), and unsteady motion of a particle in a viscous fluid, i. e. the 
Basset problem (§2). In the former analysis fractional calculus leads us to introduce 
intermediate models of viscoelasticity which generalize the classical spring-dashpot 
models. The latter analysis induces us to introduce a hydrodynamic model suitable 
to revisit in §3 the classical theory of the Brownian motion, which is a relevant 
topic in statistical mechanics. By the tools of fractional calculus we explain the long 
tails in the velocity correlation and in the displacement variance. In §4 we consider 
the fractional diffusion-wave equation, which is obtained from the classical diffusion 
equation by replacing the first-order time derivative by a fractional derivative of order 
(3 with < (3 < 2 . Led by our analysis we express the fundamental solutions (the 
Green functions) in terms of two interrelated auxiliary functions in the similarity 
variable, which turn out to be of Wright type (see Appendix), and to distinguish 
slow-diffusion processes (0 < f3 < 1) from intermediate processes (1 < f3 < 2). 

2000 Mathematics Subject Classification: 26A33, 33E12, 44A20, 45J05, 45K05, 
60E07, 60G18, 60J60, 60J65, 74D05, 76Dxx. 

This research was partially supported by the Ministry of University and by the 
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Gorenflo for fruitful discussions and comments. 



292 Fractional Calculus: Some Basic Problems 



1. LINEAR VISCOELASTICITY AND FRACTIONAL CALCULUS 

1.1 Fundamentals of Linear Viscoelasticity 

Viscoelasticity is a property possessed by bodies which, when deformed, exhibit 
both viscous and elastic behaviour through simultaneous dissipation and storage of 
mechanical energy. Here, for simplicity, we are restricting the discussion only to the 
scalar case, i.e. to one-dimensional problems. We denote the stress by a = a(x, t) and 
the strain by e = e(x, t) where x and t are the space and time variables, respectively. 

According to the linear theory of viscoelasticity, at a fixed position, the body may 
be considered a linear system with the stress (or strain) as the excitation function 
(input) and the strain (or stress) as the response function (output). Consequently, 
the response functions to an excitation expressed by the Heaviside step function 0(£) 
are known to play a fundamental role both from a mathematical and physical point 
of view, see e.g. Gross [1], Bland [2], Caputo & Mainardi [3], Christensen [4] and 
Pipkin [5]. 

We denote by J(t) the strain response to the unit step of stress (creep test), and 
by G(t) the stress response to a unit step of strain (relaxation test). These functions 
J(t) , G(t) are usually referred to as the creep compliance and relaxation modulus 
respectively, or, simply, the material functions of the viscoelastic body. In view of 
the causality requirement, both the functions are causal (i.e. vanishing for t < 0). 
The limiting values of the material functions for t — > + and t — > +oo are related 
to the instantaneous (or glass) and equilibrium behaviours of the viscoelastic body, 
respectively. As a consequence, it is usual to denote J g := J(0 + ) the glass compliance, 
J e := J(+oo) the equilibrium compliance, and G g := G(0 + ) the glass modulus, G e := 
G(+oo) the equilibrium modulus. As a matter of fact, both the material functions 
are non-negative. Furthermore, for < t < +oo , J(t) is a differentiable increasing 
function of time, i.e. 

t E 1R + , ^ > ^ < J(0+) < J(t) < J(+oo) < +oo , 

(Jib 

while G(t) is a differentiable decreasing function of time, i.e. 

dC 

t e 1R + , — < +oo > G(0+) > G(t) > G(+oo) > . 

(JLL 

The above characteristics of monotonicity of J(t) and G(t) are related respectively to 
the physical phenomena of strain creep and stress relaxation, which are experimentally 
observed. Later on, we shall outline more restrictive mathematical conditions that the 
material functions must usually satisfy to agree with the most common experimental 
observations. 
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By using the Boltzmann superposition principle, the general stress-strain relation 
can be expressed in terms of one material function [J(t) or G(t)] through a linear 
hereditary integral of Stieltjes type, namely 

e(t)= I J{t-r)da{r), or a(t) = I G(t-r)de(r). (1.1) 

J — oo J — oo 

Usually, the viscoelastic body is quiescent for all times prior to some starting instant 
that we assume as t = . Thus, under the assumption of causal histories, differentiable 
for t E M + , the representations (1.1) reduce to 

e(t)= [ J(t-T)da(r) = a(0 + )J(t)+ I J(t - r) &(r) dr , (1.2a) 
Jo- Jo 

a(t)= [ G(t - t) de(r) = e(0+) G(t) + I G(t - r) e(r) dr , (1.2b) 
Jo~ Jo 
where the superposed dot denotes time-differentiation. The lower limits of integration 

in Eqs (1.2) are written as 0~ to account for the possibility that a(t) and/or e(t) 
exhibit jump discontinuities at t = 0, and therefore their derivatives &(t) and e(t) 
involve a delta function 5(t) . Another form of the constitutive equations can be 
obtained from Eqs (1.2) integrating by parts: 

e(t) = J g a(t) + f j{t-T)a{r)dT , (1.3a) 
Jo 

a(t) = G g e(t) + [ G(t- t) e(r) dr . (1.36) 
Jo 

Here we have assumed J g > and J g < oo , see (1.7). The causal functions J(t) 
and G(t) are referred to as the rate of creep (compliance) and the rate of relaxation 
(modulus), respectively; they play the role of memory functions in the constitutive 
equations (1.3). Being of convolution type, equations (1.2) and (1.3) can be conve- 
niently treated by the technique of Laplace transforms to yield 

e(s) = sJ(s)a(s), a(s) = s G(s)e(s) . (1.4) 

Since the creep and relaxation integral formulations must agree with one another, 
there must be a one-to-one correspondence between the relaxation modulus and the 
creep compliance. The basic relation between J(t) and G(t) is found noticing the 
following reciprocity relation in the Laplace domain, deduced from Eqs (1.4), 

s J(s) = <=► J(s) G(s) = 1 . (1.5) 

sG(s) s 

Then, inverting the R.H.S. of (1.5), we obtain 

J(t) * G(t) := / J(t - t) G(t) dr = t. (1.6) 
Jo 
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Furthermore, in view of the limiting theorems for the Laplace transform we can 
deduce from the L.H.S of (1.5) that 

1 _ 1 



Ja 



Gr 



G e 



(1.7) 



with the convention that and +oo are reciprocal to each other. These remarkable 
relations allow us to classify the viscoelastic bodies according to their instantaneous 
and equilibrium responses. In fact, we easily recognize four possibilities for the lim- 
iting values of the creep compliance and relaxation modulus, as listed in Table I. 



Type 


J 9 


Je 


G g 


G e 


I 


> 


< oo 


< oo 


> 


II 


> 


= oo 


< oo 


= 


III 


= 


< oo 


= oo 


> 


IV 


= 


= oo 


= oo 


= 



Table I: The four types of viscoelasticity 
From a mathematical point of view the material functions turn out to be of the 
following form [1] 



J(t) = J g + x+J o Re(r) (l-e _ */ r ) dr + J+t, 

roc 

G{t) = G e + X - / Ra(r)e- t / T dr + G_S(t). 
Jo 



(1.8) 



where all the coefficients and functions are non negative. The function R £ (r) is 
referred to as the retardation spectrum while R a (r) as the relaxation spectrum. For 
the sake of convenience we shall denote by R*(t) anyone of the two spectra. The 
spectra must necessarily be locally summable in 1R + ; if they are summable, the 
supplementary normalization condition J Q R*(t) dr = 1 is required for the sake 
of convenience. We devote particular attention to the integral contributions to the 
material functions (1.8), i.e. 



tt(t) 



f°° / 4-1 \ d n ^> 

X+ J o Rz(t) (l-e-*/T) dr =► (-!)»_< 0, n e N . 



*(t) := X - 



R a (r) e _t / r rfr 



(-1) 



n >0, 



(1.9) 



nelN. 



, dV 
The positive functions ^(t) and $(t) are simply referred to as the creep and relaxation 
functions, respectively. According to standard definitions, see e.g. [6], the alternating 
sign properties outlined in the R.H.S. of (1.9) imply that the creep function is of 
Bernstein type, and the relaxation function is completely monotone. In particular, 
we recognize that ^(t) is an increasing function with \I/ (0) = and \I/(+oo) = x+ ° r 
+oo , while $(t) is a decreasing function with $(0) = X- or +°° an d $(+oo) = . 
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1.2 The Mechanical Models 

To get some feeling for linear viscoelastic behaviour, it is useful to con- 
sider the simpler behaviour of analog mechanical models. They are con- 
structed from linear springs and dashpots, disposed singly and in branches 
of two (in series or in parallel), as it is indicated in Fig. 1-1. 



A 



1 



i 



a) b) 



c) 



Fig. 1-1 

The elements of the mechanical models: a) Hooke, b) Newton, c) Voigt, d) Maxwell 

As analog of stress and strain, we use the total extending force and the total 
extension. We note that when two elements are combined in series [in parallel] , their 
compliances [moduli] are additive. This can be stated as a combination rule: creep 
compliances add in series, while relaxation moduli add in parallel. 

The mechanical models play an important role in the literature which is justified 
by the historical development. In fact, the early theories were established with the 
aid of these models, which are still helpful to visualise properties and laws of the 
general theory, using the combination rule. 

Now, it is worthwhile to consider the simplest mechanical models and provide 
their governing stress-strain relations along with the related material functions. We 
point out that the technique of Laplace transform allows one to easily obtain the 
requested material functions from the governing equations. 

The spring, see Fig. 1-la), is the elastic (or storage) element, as for it the force is 
proportional to the extension; it represents a perfect elastic body obeying the Hooke 
law (ideal solid). This model is thus referred to as the Hooke model. We have 

J(t) = 1/m 
G(t) = m 



a(t) = me(t) Hooke 



(1.10) 



296 Fractional Calculus: Some Basic Problems 



Newton 



The dashpot, see Fig. 1-lb), is the viscous (or dissipative) element, the force being 
proportional to the rate of extension; it represents a perfectly viscous body obeying 
the Newton law (perfect liquid). This model is thus referred to as the Newton model. 
We have 

' J(t) = t/b 
G(t) =b8(t) 

We note that the Hooke and Newton models represent the limiting cases of viscoelastic 
bodies of type / and IV, respectively. 

A branch constituted by a spring in parallel with a dashpot is known as the Voigt 
model, see Fig. 1-lc). We have 



, de 



(1.11) 



a(t) = me(t) + 6^ 



Voigt 



J(t) = - r 
m 



1 -e"*M] 
G(t) =m + b6(t) 



(1.12) 



where r e = b/m is referred to as the retardation time. 

A branch constituted by a spring in series with a dashpot is known as the Maxwell 
model, see Fig. 1-ld). We have 



a(t) + a 



da 
~dt 



de 
dt 



Maxwell 



w s a t 

j v=b + b 



G(t) 



-t/r a 



(1.13) 



where r a = a is is referred to as the the relaxation time. 

The Voigt and the Maxwell models are thus the simplest viscoelastic bodies of 
type III and II, respectively. The Voigt model exhibits an exponential (reversible) 
strain creep but no stress relaxation; it is also referred to as the retardation element. 
The Maxwell model exhibits an exponential (reversible) stress relaxation and a linear 
(non reversible) strain creep; it is also referred to as the relaxation element. 

Adding a spring either in series to a Voigt model, see Fig. l-2a), or in parallel 
to a Maxwell model, see Fig. l-2b), means, according to the combination rule, to 
add a positive constant both to the Voigt-like creep compliance and to the Maxwell- 
like relaxation modulus so that we obtain J g > and G e > . Such a model was 
introduced by Zener [7] with the denomination of Standard Linear Solid (S.L.S.). 
We have 



d 

1 + a dt 



a(t) = 



. d 
m + b — 
dt 



e(t) SLS 



J(t) = J g + x+ [l - 
G(t) =G e + X - e^/^ 



-t/r ( 



(1.14) 
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Jo 



a 
b 



X+ 



G e =m, x- = 



1 


a 


b 




~b' 




m 


m 


b 








— m , 


T a = a 


a 







(1.15) 



We point out that the condition < m < b/a ensures that x+ > X- are positive and 
hence < J g < J e < oo , < G e < G g < oo and < r a < r e < oo . The S.L.S. is the 
simplest (3-parameter) viscoelastic body of type I . On the other hand, adding a dash- 
pot either in series to a Voigt model, see Fig. l-2c), or in parallel to a Maxwell model, 
see Fig. l-2d), we obtain the simplest (3-parameter) viscoelastic body of type IV . 



1 



I 



1 



1 



I 



a) 



b) 



c) 



d) 



Fig. 1-2 

a) spring in series with Voigt, b) spring in parallel with Maxwell; 
c) dashpot in series with Voigt, d) dashpot in parallel with Maxwell. 
Based on the combination rule, we can construct models whose material functions 
are of the following type 

e -t/T e ,n 



J(t) = j g + y j n \i 



n 



+ J+t, 



G(t) = G e + J2°n e^/^" + G_ 6(t) , 



(1.16) 



where all the coefficient are non-negative. These functions must be interrelated be- 
cause of the reciprocity relation (1.5) in the Laplace domain. Appealing to the theory 
of Laplace transforms [2] , it turns out that stress-strain relation must be a linear dif- 
ferential equation with constant (positive) coefficients of the following form 



1 + l^ a «^k 



k=i 



a(t) = 



q d k 



k=i 



e(t), p = q or p = q + l. (1.17) 



Eq. (1.17) is referred to as the operator equation for the mechanical models. 
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1.3 The Fractional Viscoelastic Models 

Let us now consider a creep compliance of the form 

t a 

J<+) = #(£) = a— -, a>0, 0<a<l, (1.18) 

1 (1 + a) 

where V denotes the Gamma function. Such behaviour is found to be of some interest 
in creep experiments; usually it is referred to as power-law creep. This law appears 
compatible with the mathematical theory presented in the previous sub-section, in 
that there exists a corresponding retardation spectrum, locally summable, which reads 

fl<w = fi^4_. (1 . 19) 

7T T 

For such a model the relaxation modulus can be derived from the reciprocity relation 
(1.5) and reads 

G(t) = Q(t) = b y 6 = l/a>0. (1.20) 
1 (1 — a) 

However, the corresponding relaxation spectrum does not exist in the ordinary sense, 
in that it would be 

sinTra 1 

and thus not locally summable. The stress-strain relation in the creep representation, 
obtained from (1.1) and (1.18) is therefore 

< t ) = wn^\ f (t-r)"^- (1.22) 

r(i + a) J_ 00 

Writing da = &(t) dr and integrating by parts, we finally obtain 

e(f) = -f- f (t- r)- 1 a(r) dr = a [a(t)} , (1.23) 

where J"^ denotes the fractional integral of order a with starting point — oo , see 
Gorenflo & Mainardi [8]. 

In the relaxation representation the stress-strain relation can be obtained from 
(1.1) and (1.20). Writing de = e(r) dr , we obtain 



r( 

where 



a(t) = ^- [' («-r)- rdf(T ' 
1-a) J_ 00 



dr 



dT = b ^, (1.24) 



d a d 

±- = D™_ = Ji- Q ^ (1.25) 
dt a 00 00 

denotes the Caputo fractional derivative of order a with starting point — oo , see 
Gorenflo and Mainardi [8]. 
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Of course, for causal histories, the starting point of the integrals in (1.22-25) is 
, so that we must consider the operators J a and D™ . Since in the limit as a — > 1 
the fractional integral and derivative tend to the ordinary integral and derivative, 
respectively, we note that the classical Newton model can be recovered from (1.23) 
and (1.24) by setting a = 1 . 

In textbooks on rheology the relation (1.24), when expressed with the fractional 
derivative, is usually referred to as the Scott-Blair stress-strain law from the name of 
the scientist [9] , who in earlier times proposed such a constitutive equation to intro- 
duce a material property that is intermediate between the elastic modulus (Hooke 
solid) and the coefficient of viscosity (Newton fluid). 

The use of fractional calculus in linear viscoelasticity leads to a generalization 
of the classical mechanical models in that the basic Newton element (dashpot) is 
substituted by the more general Scott-Blair element. In fact, we can construct the 
class of these generalized models from Hooke and Scott-Blair elements, disposed 
singly and in branches of two (in series or in parallel). The material functions are 
obtained using the combination rule; their determination is made easy if we take into 
account the following correspondence principle between the classical and fractional 
mechanical models, as stated by Caputo & Mainardi [3], 

t a 



t ->■ 



(0 < a < 1) < 



r(i + a) ' 



(1.26) 



r(i-a) ' 

e -0-_> Ea [-(f/T) a ], 

where E a denotes the Mittag-Lemer function of order a , heavily used in [8] . 

We verify the correspondence principle by considering the fractional S.L.S., for- 
merly introduced by Caputo & Mainardi [10] in 1971. Such model is based on the 
following operator equation of fractional order, which generalises the operator equa- 
tion (1.14) for the S.L.S., 



1 + a 



d' 



dt a 



a(t) = 



m + b 



d l 



dt 1 



c(t) , < a < 1 . 



This equation is better analysed in the Laplace domain where we obtain 

1 l + as a 



(l + as a ) a(s) = {m + bs a ) e(s) 



sJ(s) = 



(1.27) 



(1.28) 



sG(s) m + bs a 

From the fractional operator equation we can obtain as particular cases, besides the 
trivial elastic model (a = b = 0) and the fractional Newton or Scott-Blair model 
(a = m = , b = fH) already considered, the fractional Voigt model (a = 0) and the 
fractional Maxwell model (m = 0). 
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Working in the Laplace domain and then inverting, we obtain for the fractional 
Voigt and Maxwell models 



a(t) = me(t) + b 



d a e 
dt* 



Fractional Voigt 



J(t) = I{l-E Q [-(t/r e ) a ]} 

fit 

+—0, 

G(t) = m + b 



(1.29) 



r(i -a) 



. . d a a , d a e , 
o(t) + a —r— = b t— r r actional Maxwell 



, « Oj 1 t a 
J(t) = T + T 



dt a 



dt c 



b bT(l + a) 

G(t) = -E a [-{t/r a r\ 
a 



(1.30) 



where (r e ) a = b/m and (T a ) a = a. 

Having recognized with (1.29-30) the validity of the Caputo-Mainardi correspon- 
dence principle for the basic models, we are allowed to use this principle to obtain the 
material functions of higher models, including the fractional S.L.S., along with the 
corresponding operator equations of fractional order. Thus, by generalizing (1.16), 
we obtain 

J(t) = J g + J2Jn{l-E a [-(f/r eiB ) a ]} + J+ *° , 

t -a (1-31) 

G(t) = G e + J2G n E a [-(t/r a , n r] + G. — — — , 

n ^ ' 

where all the coefficients are non negative. Extending the procedures of the classical 
mechanical models, we will get the fractional operator equation in the form which 
properly generalises (1.17), i.e. 



fc=i 



d ak 
dt a K 



a(t) 



rn 



d 



Oik 



k=l 



dt a " 



e(t), a k = k + a-l. (1.32) 



We conclude this section pointing out the presence of the Mittag-Leffler function 
in (1.31). In fact, the creep and relaxation functions for the fractional models contain 
contributions of type 



Re(r) (l-e^ /T ) dr, 

m = X-E a [-(t/r a r]=X- / RMe-^dr. 

Jo 



(1.33) 



Denoting as usual by * the suffix e or a, the analytical expressions of the retardation 
and relaxation spectra turn out to be identical, namely 



R*(t) = — 



sin an 



tit (r/r*) a + (r/r*) a + 2 cos aix 



(1.34) 
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This result can be deduced from the spectral representation of the Mittag-Leffler func- 
tion E a [— (t/r*)"], as shown by Caputo and Mainardi [3], and recently by Gorenflo & 
Mainardi [8] in the framework of their analysis of the fractional relaxation equation. 

We can have a better insight of the spectral function R*(t) and of the relaxation 
function E a [— (t/r*) a ] by showing the corresponding plots for a few values of a. 
Assuming r* = 1 , we could simply refer to the plots reported in [8] by Fig. la and 
Fig. 2a, but, for the sake of convenience, we prefer to exhibit them again in Fig. 1-3 
and Fig. 1-4, hereafter. 



2 



1.5 



0.5 



1 R.(x) 








a=0.90 






^cx-0.75 >y 






^^a=0.50 




a=CL25 





0.5 1 1.5 2 



Fig. 1-3 

Spectral function i2„(r) for a = 0.25 , 0.50 , 0.75 , 0.90 . 

From the plots of R*(r) in Fig. 1-3 we can easily recognize the effect of the 
variation of a on the character of the spectral function; for a — > 1 the spectrum 
becomes sharper and sharper until for a = 1 it reduces to be discrete with a single 
retardation/relaxation time. We also recognize that i?*(r) is a decreasing function 
of r for < a < a* where a* ~ 0.736 is the solution of the equation a = sin an ; 
subsequently, with increasing a , it first exhibits a minimum and then a maximum 
before tending to the impulsive function S(r — r*) as a — > 1 . Recalling the analysis of 
the fractional relaxation equation by Gorenflo and Mainardi [8], we recognize that, 
compared to the exponential obtained for a = 1 , the fractional relaxation function 
exhibits very different behaviours, as can be seen from the plots of E a (—t a ) in Fig. 
1-4. In particular, we point out the leading asymptotic behaviours at small and large 
times, 

v ' \t a /r(l-a), as t-^+oo. v ' 

Compared to the solution exp(— t) for the classical models (a = 1), the solution 
E a (— t a ) for the fractional models (0 < a < 1) exhibits initially a much faster decay 
(the derivative tends to — oo in comparison with —1), and for large times a much 
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slower decay (algebraic decay in comparison with exponential decay) . In view of its 
final slow decay, the phenomenon of fractional relaxation is usually referred to as a 
super-slow process. 



1 E (-t a ) 








_ ot-1 


a=0.25 

ot=0.50 

a=0.75 



5 10 t 15 



Fig. 1-4 

Relaxation function E a (-t a ) for a = 0.25 , 0.50 , 0.75 , 1 . 

1.4 Bibliographical remarks 

A number of authors have, implicitly or explicitly, used fractional calculus as an 
empirical method of describing the properties of viscoelastic materials. 

In the first half of this century Gemant [11-12] and, later, Scott-Blair [9, 13] 
were early contributors in the use of fractional calculus to study phenomenological 
constitutive equations for viscoelastic media. 

Independently, in the former Soviet Union, Rabotnov [14-15] introduced his the- 
ory of hereditary solid mechanics with weakly singular kernels, that implicitly requires 
fractional derivatives. This theory was developed also by other soviet scientists in- 
cluding Meshkov and Rossikhin, see e.g. [16], and Lokshin and Suvorova, see e.g. 
[17]. 

In 1971, extending earlier work by Caputo [18-20], Caputo and Mainardi [3,10] 
suggested that derivatives of fractional order could be successfully used to model the 
dissipation in seismology and in metallurgy. Since then up to nowadays, applications 
of fractional calculus in rheology have been considered by several authors. Without 
claim of being exhaustive, we now quote some papers of which the author became 
aware during the last 25 years. In addition to Caputo [21-24] and Mainardi [25- 
26] we like to refer to Smith and de Vries [27], Scarpi [28], Stiassnie [29], Bagley 
and Torvik [30-33], Rogers [34], Koeller [35-36], Koh and Kelly [37], Friedrich [38], 
Nonnenmacher and Glockle [39-40], Makris and Constantinou [41], Heymans and 
Bauwens [42], Schiessel & al [43], Gaul & al [44], Beyer and Kempfle [45], Fenander 
[46], Pritz [47], Rossikhin & al [48-49], and Lion [50]. 



F. Mainardi 303 



2. THE BASSET PROBLEM VIA FRACTIONAL CALCULUS 
2.1 Introduction 

The dynamics of a sphere immersed in an incompressible viscous fluid represents a 
classical problem, which has many applications in flows of geophysical and engineering 
interest. Usually, the low Reynolds number limit (slow motion approximation) is 
assumed so that the Navier-Stokes equations describing the fluid motion may be 
linearised. 

The particular but relevant situation of a sphere subjected to gravity was first 
considered independently by Boussinesq [51] in 1885 and by Basset [52] in 1888, who 
introduced a special hydrodynamic force, related to the history of the relative accel- 
eration of the sphere, which is nowadays referred to as Basset force. The relevance 
of these studies was in that, up to then, only steady motions or small oscillations 
of bodies in a viscous liquid had been considered starting from Stokes' celebrated 
memoir on pendulums [53], in 1851. The subject matter was considered with more 
details in 1907 by Picciati [54] and Boggio [55], in some notes presented by the great 
Italian scientist Levi-Civita. The whole was summarised by Basset himself in a later 
paper [56], and, in more recent times, by Hughes and Gilliand [57]. 

Nowadays the dynamics of impurities in unsteady flows is quite relevant as shown 
by several publications, whose aim is to provide more general expressions for the 
hydrodynamic forces, including the Basset force, in order to fit experimental data 
and numerical simulations, see e.g. [58-66]. 

In the next section we shall recall the general equation of motion for a spherical 
particle, in a viscous fluid, pointing out the different force contributions due to effects 
of inertia, viscous drag and buoyancy. In particular, the so-called Basset force will 
be interpreted in terms of a fractional derivative of order 1 /2 of the particle velocity 
relative to the fluid. Based on our recent works [67-68], we shall introduce the 
generalized Basset force, which is expressed in terms of a fractional derivative of 
any order a ranging in the interval < a < 1 . This generalization, suggested by a 
mathematical speculation, is expected to provide a phenomenological insight for the 
experimental data. 

In section §2.3 we shall consider the simplified problem, originally investigated 
by Basset, where the fluid is quiescent and the particle moves under the action of 
gravity, starting at t = with a certain vertical velocity. For the sake of generality, 
we prefer to consider the problem with the generalized Basset force and will provide 
the solution for the particle velocity in terms of Mittag-Leffler -type functions. The 
most evident effect of this generalization will be to modify the long-time behaviour 
of the solution, changing its algebraic decay from t -1 / 2 to t~ a . This effect can be of 
some interest for a better fit of experimental data. 
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2.2 The Equation of Motion 

Let us consider a small rigid sphere of radius tq, mass m p , density p p , initially 
centred in X(t) and moving with velocity V(t) in a homogeneous fluid, of density 
Pf and kinematic viscosity v , characterized by a flow field u(x,t) . In general the 
equation of motion is required to take into account effects due to inertia, viscous 
drag and buoyancy, so it can be written as 

m p ^ = F l + F d + F g , (2.1) 

where the forces on the R.H.S. correspond in turn to the above effects. According to 
Maxey and Riley [60] these forces read, adopting our notation, 



Fi = m f 



Du 
~Dt 



1 



TO/ 



F d = — 
p 

where to/ = 
particle, and 



X(t) 

[V(t)-u(X(t),t)} + 





Du 

Dt X (t) 
d[V(r)-u(X(r),r)]/dr 



dr 



F g = (mp -m f )g. 



(2.2) 

(2.3) 
(2.4) 



(4/3)7rroP/ denotes the mass of the fluid displaced by the spherical 



to := — , 
v 



1 a 9 -i 
-:=Q 7V r iypf = -m f T . 



(2.5) 
(2.6) 



The time constant tq represents a sort of time scale induced by viscosity, whereas 
the constant p is usually referred to as the mobility coefficient. 

In (2.2) we note two different time derivatives, D/Dt , d/dt , which represent the 
time derivatives following a fluid element and the moving sphere, respectively, so 



Du 
TJt 



x(t) 



du 



+ (it • V) u(x,t) 



d 
dt 



u[X(t),t] 



du 



+ (V-V)u(x,t) 



where the brackets are computed at x = X(t) . 

The terms on the R.H.S. of (2.2) correspond in turn to the effects of pressure 
gradient of the undisturbed flow and of added mass, whereas those of (2.3) represent 
respectively the well-known viscous Stokes drag, that we shall denote by Fs , and to 
the augmented viscous Basset drag denoted by Fb ■ Using the characteristic time To, 
the Stokes and Basset forces read respectively 



F s = --m f r - 1 [V(t)-u(X(t),t)], 



9 -1/2 
o m / T 



1 



d[V(T)-u(X(T),T)]/dT 



dr 



(2.7) 
(2.8) 
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We thus recognize that the time constant tq provides the natural time scale for 
the diffusive processes related to the fluid viscosity, and that the integral expression 
in brackets at the R.H.S. of (2.8) just represents the Caputo fractional derivative of 
order 1/2 , with starting point — oo , of the particle velocity relative to the fluid *. 

We now introduce the generalized Basset force by the definition 

n = -\rn f rr 1 ^[V{t)-u{X{tlt)], < a < 1 , (2.9) 

where the fractional derivative of order a is in Caputo's sense, in agreement with the 
notation introduced in §1.3 for the fractional viscoelastic models, see (1.25). 
Introducing the so-called effective mass 

1 

m e := m p + -m f , (2.10) 

and allowing for the generalized Basset force in (2.3), we can re- write the equation of 
motion (2.1-4) in the more compact and significant form, 



dV 3 Du 9 
me ^t = 2 mf Dt ~ 2 mf 



1 1 d a 
+ 



T{) ^ TQ i- a dfa ^ (V-u) + (m p -m f )g, (2.11) 

that we refer to as the generalized equation of motion. Of course, if in (2.11) we put 
a = 1/2 , we recover the basic equation of motion with the original Basset force. 

2.3 The (Generalized) Basset Problem 

Let us now assume that the fluid is quiescent, namely u(x, t) = , Vx, t , and the 
the particle starts to move under the action of gravity, from a given instant to = 
with a certain velocity V(0 + ) = Vq , in the vertical direction. This was the problem 
considered by Basset [52] , that was first solved by Boggio [55] , in a cumbersome way, 
in terms of Gauss and Fresnel integrals. 

Introducing the non-dimensional quantities (related to the densities pf , p p of the 
fluid and particle), 

X'-=^, /3:=^— = -?-, (2.12) 
Pf 2p P + Pf 1 + 2% 1 ' 

we find it convenient to define a new characteristic time 

a e := [im e = r //3, (2-13) 
see (2.5), (2.10), (2.12), and a characteristic velocity (related to the gravity), 

V s = (2/9)(x-l)gr Q . (2.14) 



* Presumably, the first scientist who has pointed out the relationship between the Basset force 
and the fractional calculus has been Tatom [69] in f 988. However, Tatom has limited himself to note 
this fact, without treating any related problem by the methods of fractional calculus. 
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1 + T0 d^ 



V + — V s . (2.15) 



Then we can eliminate the mass factors and the gravity acceleration in (2.11) and 
obtain the equation of motion in the form 

dV _ 1_ 

dt a e 

If the Basset term were absent, we obtain the classical Stokes solution 

V(t) = V S + (Vb - V s ) e-t/** , (2.16) 

where a e represents the characteristic time of the motion, and Vs the final value 
assumed by the velocity. Later we shall show that in the presence of the Basset term 
the same final value is still attained by the solution V(t), but with an algebraic rate, 
which is much slower than the exponential one found in (2.16). 

In order to investigate the effect of the (generalized) Basset term, we compare 
the exact solution of (2.15) with the Stokes solution (2.16); for this aim we find 
it convenient to scale times and velocities in (2.15) with {a e , Vs}, i.e. to refer to 
the non dimensional quantities t' = t/a e , V = V/Vs , V ' = Vq/Vs ■ The resulting 
equation of motion reads (suppressing the apices) 



d d a 

-r + a- hi 

dt dt a 



V(t) = l, V(0 + ) = V , a = /T>0, 0<a<l. (2.17) 



This is the composite fractional relaxation equation treated by Gorenflo and Mainardi 
[8] in §4.1 by using the Laplace transform method. Recalling that in an obvious 
notation we have 

V(t) - V(s) , ^V(t) - s a V(s) - s*- 1 V , < a < 1 , (2.18) 
the transformed solution of (2.17) reads 

V{s) = M(s)V +-N(s), (2.19) 
s 



where 



1 _|_ n s a-l „ 1 

M(s) = — -, N{s) = -. (2.20) 

' s + as a + l w s + as a + l v ' 



Noting that 

- N(s) = - - Mis) -r / N(r)dr= 1-M(t) <^ N(t) = -M'(t), (2.21) 
s s Jo 

the actual solution of (2.17) turns out to be 

V(t) = l + (V -l)M(t), (2.22) 

which is "similar" to the Stokes solution (2.16) if we consider the substitution of e~* 
with the function M(t) . 
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In [67-68] Mainardi, Pironi and Tampieri have used a factorisation method to 
invert N(s) and henceforth M(s) , using a procedure indicated by Miller and Ross 
[69], which is valid when a is a rational number, say a = p/q, where p, q G IN , p < q . 
In this way the actual solution can be finally expressed as a linear combination of 
certain incomplete gamma functions. This algebraic method is of course convenient 
for the ordinary Basset problem (a = 1/2), but becomes cumbersome for q > 2 . 

Here, following the analysis in [8], we prefer to adopt the general method of 
inversion based on the complex Bromwich formula. By this way we are free from the 
restriction of being a a rational number and, furthermore, we are able to provide an 
integral representation of the solution, convenient for numerical computation, which 
allows us to recognize the monotonicity properties of the solution without need of 
plotting. 

We now resume the relevant results from [8] using the present notation. The 
integral representation for M(t) turns out to be 

poo 

M(t) = / e~ rt K(r)dr, (2.23) 
Jo 

where 

* . 1 ar a_1 sin (mr) 

K ( r ) = ~ T, T5 5-5 TTFr T 1 V > • 2 - 24 

n (1- r) 2 + a 2 r 2a + 2(1- r) a r a cos (an) 

Thus M(i) is a completely monotone function [with spectrum K(r)], which is de- 
creasing from 1 towards as t runs from to oo . The behaviour of M(t) as t — > 0+ 
and t — ?■ oo can be inspected by means of a proper asymptotic analysis, as follows. 

The behaviour as t — > + can be determined from the behaviour of the Laplace 
transform M(s) = s -1 - s~ 2 + O (s~ 3+a ) , as Re {s} -> +oo . We obtain 

M(t) = 1 -t + (t 2 " a ) , as t^0+. (2.25) 

The spectral representation (2.23-24) is suitable to obtain the asymptotic be- 
haviour of M(t) as f -> +00 , by using the Watson lemma. In fact, expanding 
the spectrum K(r) for small r and taking the dominant term in the corresponding 
asymptotic series, we obtain 

M(t)~a /"° . =a sin ( a7r ) [°° e-rt r «-i dr ^ as t ^ ^ . (2 . 26) 
r(l-a) 7T J Q 

Furthermore, we recognize that 1 > M(t) > > , < t < oo , namely, the 
decreasing plot of M(t) remains above that of the exponential, as t runs from to 
oo . Although both the two functions tend monotonically to , the difference between 
the two plots increases with t: at the initial point t = , both the curves assume 
the unitary value and decrease with the same initial rate, but as t — > oo they exhibit 
very different decays, algebraic (slow) against exponential (fast). 
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For the ordinary Basset problem it is convenient to report the result obtained 
by the factorisation method [67-68]. In this case we must note that a = y/]3 , see 
(2.17), ranges from to 3 since from (2.12) we recognize that /3 runs from (x = oo , 
infinitely heavy particle) to 9 (x = , infinitely light particle) . 

The actual solution is obtained expanding M(s) into partial fractions and then 
inverting. Considering the two roots X± of the polynomial P(z) = z 2 + a z + 1 , with 
z = s 1 / 2 we must treat separately the following two cases 

z) < a < 2 , or 2<a<3, and ii) a = 2 , 

which correspond to two distinct roots (A+ 7^ A_), or two coincident roots (A+ = 
A_ = —1), respectively. We obtain 

[S) ~ s + asW + l ~ S V2 ( s i/2 _ A+ ) + s i/2( s i/2_ A_) ' ( ] 



with 



U) a = 2 ^ 13 = 4, x = 5/8, 

~ 1 + 2 s -1 / 2 1 2 

M ^ = 8 + 2 s l/2 + 1 = ( s l/2 + !)2 + s l/2 ( s l/2 + 1)2 ' ( 2 " 29 ) 

The Laplace inversion of (2.27 — 29) can be expressed in terms of Mittag-Leffler 
functions of order 1/2 , Ei/ 2 (\\/i) = exp(A 2 t) erfc(— Ay^) , as shown in the Appendix 
of [8]. We obtain 

{ ii) (1 - 2t) E U2 (-Vt) + 2 v^V . 

We recall that the analytical solution to the classical Basset problem was formerly 
provided by Boggio [55] in 1907 with a different (cumbersome) method. One can show 
that our solution (2.30), derived by the tools of the Laplace transform and fractional 
calculus, coincides with Boggio's solution. Also Boggio arrived at the analysis of the 
two roots A± but his expression of the solution in the case of two conjugate complex 
roots (x > 5/8) given as a sum of Fresnel integrals could induce one to forecast 
unphysical oscillations, in the absence of numerical tables or plots. This disturbed 
Basset who, when he summarised the state of art about his problem in a later paper 
of 1910 [56], thought there was some physical deficiency in his own theory. With 
our integral representation of the solution, see (2.23-24), we can prove the monotone 
character of the solution, even if the arguments of the exponential and error functions 
are complex. 
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In order to have some insight about the effects of the two parameters a and a on 
the (generalized) Basset problem we exhibit some (normalized) plots for the particle 
velocity V(t), corresponding to the solution of Eq. (2.17), assuming for simplicity a 
vanishing initial velocity (Vo = 0). 

We consider 3 cases for a, namely a = 1/2 (the ordinary Basset problem) and 
a = 1/4, 3/4 (the generalized Basset problem), corresponding to Figs 2-1, 2-2, 2-3, 
respectively. For each a we consider four values of a corresponding to x '■= Pp/Pf = 
0.5, 2, 10, 100 . For each couple {a , x} we compare the Basset solution (in continuous 
line) with its asymptotic expression (in dashed-dotted line) for large times and the 
Stokes solution (dashed line). We remind that the Stokes solution is the solution of 
Eq. (2.17) with a = and hence is independent of a . 

From these figures we can recognize the retarding effect of the (generalized) Basset 
force, which is more relevant for lighter particles, in reaching the final value of the 
velocity. This effect is of course due to the algebraic decay of the function M(t), see 
(2.26), which is much slower than the exponential decay of the Stokes solution. 
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Fig. 2-1 

The normalized velocity V(t) for a = 1/2 and x = 0.5 , 2 , 10 , 100 
Basset exact - - ; Basset asymptotic — Stokes . 
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The normalized velocity V(t) for a = 1/4 and x = 0-5 , 2 , 10 , 100 
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The normalized velocity V(t) for a = 3/4 and x = 0.5 , 2 , 10 , 100 
Basset exact - - ; Basset asymptotic — Stokes . 
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3. BROWNIAN MOTION AND FRACTIONAL CALCULUS 
3.1 Introduction 

According to the classical approach started by Langevin normal diffusion and 
Brownian motion are associated with the Langevin equation. More specifically, the 
classical Langevin equation addresses the dynamics of a Brownian particle through 
Newton's law by incorporating the effect of the Stokes fluid friction and that of 
thermal fluctuations in the vicinity of the particle into a random force, see e.g. Wax 
[70], Fox and Uhlenbeck [71], Fox [72], Kubo et al [73]. 

Since the pioneering computer experiments by Alder and Wainwright [74] in 1970, 
which have shown that the velocity autocorrelation function for a Brownian particle 
in a dense fluid goes asymptotically as t -3 / 2 instead of exponentially as predicted 
by stochastic theory, many attempts have been made to reproduce this result by 
purely theoretical arguments, see e.g. [75-97]; in most cases hydrodynamic models 
are adopted. 

Recently, a great interest on the subject matter has been raised because of the 
possible connection among long-time correlation effects, fractional Brownian motion 
and anomalous diffusion, see e.g. [98-102]. We recall that anomalous diffusion is 
the phenomenon, usually met in disordered or fractal media, according to which the 
displacement variance is no longer linear in time but proportional to a power a of 
time with < a < 1 (slow diffusion) or 1 < a < 2 (fast diffusion), see Bouchaud and 
Georges [99] for a review. 

We also point out that, in view of the linear-response theory, Kubo in 1966 [103] 
stated a, fluctuation- dissipation theorem * by introducing a generalized Langevin equa- 
tion (GLE), with an indefinite memory function as an integral kernel. In other words, 
this theorem may be represented by a stochastic equation describing the fluctuation, 
which is a generalization of the classical Langevin equation; in the GLE the friction 
force becomes retarded or frequency dependent and the random force is no longer 
a white noise. As a matter of fact, the hydrodynamic models introduced in the 
literature appear as particular cases of Kubo's GLE. 

Here, after resuming in §3.2 the classical results derived from the ordinary 
Langevin equation, in §3.3 we shall revisit a hydrodynamic model which takes into 
account, in addition to the Stokes viscous drag, the inertial effect due to the added 
mass and the retarding effect due to the Basset memory force. So doing, we obtain 
a stochastic differential equation which contains a time derivative of order 1/2. This 
GLE will be referred to as the fractional Langevin equation. 



* For a critical analysis of Kubo's fluctuation-dissipation theorem see Felderhof [104] 
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The present approach is based on a recent analysis carried out by the author and 
collaborators [105-106], in order to model the Brownian motion more realistically 
than in the classical approach (based on the Langevin equation). 

Using Kubo's fluctuation- dissipation theorem and the techniques of fractional cal- 
culus, we shall provide the analytical expressions of the autocorrelation functions 
(both for the random force and the particle velocity) and of the displacement vari- 
ance. Consequently, the well-known results of the classical theory of the Brownian 
motion will be properly generalized. 

In the final section, §3.4, we shall present and discuss some numerical results 
implied by our analysis. 

3.2 The Classical Approach to the Brownian Motion 

We assume that the Brownian particle of mass m p executes a random motion in 
one dimension with velocity V = V(t) and displacement X = X(t). The classical 
approach to the Brownian motion is based on the following stochastic differential 
equation {Langevin equation) 

mp ^- = F(t) + R(t), (3.1) 

where F(t) denotes the frictional force exerted from the fluid on the particle and R(t) 
denotes the random force arising from rapid thermal fluctuations, subjected to the 
condition ( R(t) ) = . As usual, we have denoted with brackets the average taken 
over an ensemble in thermal equilibrium. Therefore the total force has been divided 
into a mean force F and a fluctuating force R. The fact that F(t) is independent of 
the fluid variables is due to the boundary condition that the fluid velocity be equal 
to the particle velocity, V(t) , at the surface of the particle. 

Assuming for the mean force the familiar Stokes approximation for a drag of 
spherical particle of radius ro , we obtain the classical formula 

F = --V(t), -=67rr p f u, (3.2) 
p p 

where p denotes the mobility coefficient and pf and v are the density and the kine- 
matic viscosity of the fluid, respectively. In this approximation the time derivative of 
the fluid velocity field has been neglected. If we introduce the friction characteristic 
time a p := pm p , the Langevin equation (3.1) explicitly reads 



^ = --V(t) + ^-R(t) 
dt a p m p 



(3.3) 
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The stochastic processes V(t) and R(t) are assumed to be Gaussian-Markovian 
and stationary. The stationarity implies that the autocorrelation functions Cy and 
Cr depend only on the time shift, namely 

Cvitut) := (V(t 1 )V(t 1 +t)} = C v (t), (3.4) 

C R {t u t) := ( R{t x ) R(h +t)) = C R (t) , (3.5) 

for any t\ and t . Hereafter we assume t\ = and t > . 

Following the classical approach to the Brownian motion, we require that the 
variance of the velocity at t = , CV(0) = (V a2 (0) ) , satisfies the equipartition law 
for the energy distribution, i.e. 

m p (V 2 (0)} = kT <^ a p (V 2 (0) } = k T , (3.6) 

where is the Boltzmann constant, as if the Brownian particle were kept for a 
sufficiently long time in the fluid at (absolute) temperature T , and that the random 
force is uncorrelated to the particle velocity at t = , i.e. 

(V(0)R(t)) = 0, t>0. (3.7) 

As well known, the previous assumptions lead to the relevant results, 

C v (t) = (V 2 (0))e- t / a P , t>0, (3.8) 

9 

777, 

C R (t) = ^(V 2 (0))5(t), t>0, (3.9) 

where S(t) denotes the Dirac distribution. The result (3.8) shows that the velocity 
autocorrelation function decays exponentially with characteristic time a v , whereas 
(3.9) means that R(t) is a white noise. 

It can be readily shown that the mean squared displacement of the Brownian 
particle (starting at the origin at to = 0), i.e. the displacement variance, is given by 

(X 2 (t))=2 [(t-r)C v {r)dT = 2 [ dr x I Cy(r)rfr, £>0. (3.10) 
Jo Jo Jo 

For this it is sufficient to recall that X(t) = Jq V(t') dt' , and to use the definition 
(3.4) of Cy(t) for t > to = . As a consequence of (3.8) and (3.10) we obtain 

(X 2 (t)}=2a p (V 2 (0)} [t-o-p^-e-*/^)] , £>0, (3.11) 

from which we recognize that for sufficiently large times the variance increases linearly 
with time. 
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It is usual to introduce the diffusion coefficient as 

D:=lim^^. (3.12) 

t— S-oo 2 1 

Then from (3.11-12) we obtain the chain of equalities 

/•OO 

V = a p (V 2 (0)) = C v (t)dt, (3.13) 
Jo 

and, using (3.7), 

V = /jkT. (3.14) 

The identity (3.14) is known as Einstein relation. In particular, we point out the 
asymptotic behaviour of the variance for large times, 

(X 2 (t)} = 2Vt [1 - (t/o-p)- 1 +EST] , as t ->■ oo , (3.15) 

where EST denote exponentially small terms. 

3.3 The Hydrodynamic Approach to the Brownian Motion 

On the basis of hydrodynamics, the Langevin equation (3.3) is not completely 
correct, since it ignores the effects of the added mass and Basset history force, which 
are due to the acceleration of the particle. This was formerly pointed out in the early 
seventies by a number of authors, just after the cited computer experiments by Alder 
and Wainwright [74]. 

The added mass effect requires to substitute the mass of the particle with the 
so-called effective mass, m e introduced in (2.10). As a consequence, in order to keep 
unmodified the mobility coefficient in the Stokes drag, we have to introduce a new 
friction characteristic time, a e , such that 

»:=^ = ^- <=► a e :=a + with X := Et . ( 3 .16) 

m p m e V 2 X/ Pf 

The corresponding Langevin equation is obtained form (3.3) by replacing m p with 
m e and a p with a e . With respect to the classical analysis, it turns out that the 
added mass effect, if it were present alone, would be only to lengthen the time scale 
(a e > a p ) in the exponentials entering the basic formulas (3.8) and (3.11) and to 
decrease the velocity variance (V 2 (0)} , consistently with the energy equipartition law 
at the same temperature, 

m e (V 2 (0)} = kT <^ a e {V 2 (0)} =/ikT. (3.17) 



Consequently, the diffusion coefficient turns out to be not altered by the added mass 
effect and the Einstein relation still holds. 
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In view of Kubo's fluctuation- dissipation theorem, an arbitrary retarding effect 
in the friction force (in particular that due to the Basset force) can be taken into 
account by introducing a suitable memory function ^(t) in the Langevin equation. 
The consequent GLE reads (in our notation) 



dV_ 
~dt 



1 (t - t)V(t) dr + — R(t) 

m e 



t>0 



(3.18) 



where, as usual, the limits of integration are extended to account for the possibility 
of Dirac-type distributions. The fluctuation- dissipation theorem can be readily ex- 
pressed by the Laplace transforms, see e.g. Mainardi and Pironi [105]. In our notation 
this theorem leads to 

(V 2 (0)) 



C v (s) :=(V(0)V(t) 



s + 7(s) ' 



(3.19) 



and 



C R (s) : = ( R(0) R(t) ) = ml( V 2 (0) ) 7(a) . (3.20) 

The classical results are easily recovered for t > by noting that, in the absence of 
added mass and retarding effects, we get 7(5) = l/a p 4- ^(t) = 8{t)/a p . 

Taking into account both the added mass and the Basset history force (whose 
expression has been given in the previous section in terms of a fractional derivative) 
the Langevin equation (3.3) turns out to be modified into 



dV_ 
~dt 



1 



1 + V 7 ^) 



d 1 ' 2 
dt 1 / 2 



V(t) + —R(t), r : = 



m e 



(3.21) 



Here the fractional derivative is intended in the Caputo sense with starting point 
to = , i.e. 

dt^ V{t) ^ J 



dr. 



(3.22) 



Vt^f 

We agree to refer to (3.21) as the fractional Langevin equation. 

We easily recognize that our fractional Langevin equation (3.21) can be considered 
a particular case of the GLE (3.18) by noting that 



7 (s) 



1 



1 



1 



8{t)-^r -^t-^ 2 m 



(3.23) 



where 0(£) is the Heaviside step function. Therefore the expression for 7(2) turns out 
to be defined only in the sense of distributions. Specifically, S(t) is the well-known 
Dirac delta function and t~ 3 / 2 0(£) is the linear functional over test functions, <f>(t) , 
such that 

m) - 0(0)] dt _ 



For more details on distributions, see e.g. [107] or [108]. 
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The significant change with respect to the classical case results from the t~ 3 / 2 
term. Not only does it imply a non-instantaneous relationship between the force and 
the velocity, but also it is a slowly decreasing function so that the force is effectively 
related to the velocity over a large time interval. The representation of the force in 
terms of distributions, as required by the GLE, is not strictly necessary since we can 
use the equivalent fractional form. 

Let us consider the autocorrelation for the random force. The inversion of the 
Laplace transform Cr(s) yields, by (3.22-23), 



2 

mi 



c*(f) = -M^(o)> 



t>0, (3.24) 



to be compared with the classical result (3.9). Thus, we recognize that, in the presence 
of the Basset history force, the random force can no longer be represented uniquely 
by a white noise; an additional " fractional" or " coloured" noise is present due to the 
term t -3 / 2 which, as already noted, is to be interpreted in the sense of distributions. 
Since the fluctuating force is no longer uncorrelated at different times, the fractional 
Langevin equation does not represent a Markovian process. Nevertheless, it is still 
Gaussian (since the Gaussian nature of the driving sources for the fluid is assumed), 
and stationary (in view of the time-shift invariance). since 

Let us now consider the autocorrelation for the velocity field. Inserting (3.23) in 
(3.21), it turns out as 

V{ ) ~ s+[l + 0**1/3] /a e ~ s + vW^V2 + 1/ae ' {3 - 2b) 

where (3 := To/a e , see (2.12-13). We first note that the effect of the Basset force is 
expected to be negligible for f3 — > (x '•= Pp/Pf — > *- e - f° r particles which are 
sufficiently heavy with respect to the fluid. In this case we can assume also a e ~ a p 
so the classical results (3.8), (3.9) and (3.11) turn out to be true. 

A first result concerning the asymptotic behaviour of Cy (t) as t — > oo can be easily 
obtained from (3.25) by applying the asymptotic theorem for the Laplace transform 
as s — > , see e.g. Doetsch [109]. In fact, from 

C v (s) ~ a e ( V 2 (0) ) (1 - y^o~ s 1 ' 2 ) , s -+ , 

we get 

C v (t) ~ ( V 2 (0) ) vW(4^I (t/or e y 3/2 , t^oo. (3.26) 

The presence of such a long-time tail is thus in agreement with that formerly observed 
in computer simulations by Alder and Wainwright [74]. 
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The explicit inversion of the Laplace transform in (3.25) can be carried out in a 
way similar to that used in the (deterministic) ordinary Basset problem treated in 
the previous Section, see (2.27-30). For this purpose we need to consider the function 



and recognize that 

WW) = NW °> } * " S{ " S)= S + Vg^U + ■ (3 ' 28> 

Thus, the actual solution is obtained by expanding N(s) into partial fractions and 
then inverting. 

We first obtain 

o 4± + d= 

» } S V2 (s l/2_ A+) + s l/2 (s l/2_ A _ ) ' 

N{S) = s + asW + 1 = \ . , 1 (3 - 29) 

U) + l) 2 ' 

where A± and A± are given by (2.28), and the distinction of cases i) and ii) is the 
same as there. 

Then, the Laplace inversion of (3.29) can be expressed in terms of Mittag-Leffler 
functions of order 1/2 , Ei/ 2 (\Vi) = exp(A 2 t) erfc(— Ay^) , as shown in the Appendix 
of [8]. We obtain, using (3.28), 



Cy(t) 



i) A + E 1/2 (X + ^t/7 e ) + A_ E l/2 (X_ ^t/7 e ) , 

(3.30) 

ii) (l + 2t/a e )E 1/2 (-^/t/^~ e )-(2/V^)Vt/^~ e . 



Furthermore, it can be shown that N(t) is a completely monotone function for 
t > , decreasing from 1 to , as t runs from to oo . 

Let us now consider the displacement variance, which is provided by the repeated 
integral of the velocity autocorrelation as indicated in (3.10). From the Laplace trans- 
form (X 2 (s) ) = 2Cy(s)/s 2 , we first derive the asymptotic behaviour of (X 2 (t) ) as 
t — > oo . We easily obtain 

(X 2 (t)) = 2Dt{l-2 v ^(t/a e )- 1 / 2 + [(t/ae)" 1 ]} , t -)• oo , (3.31) 

where D is the diffusion coefficient defined in (3.12-14). 
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The explicit expression of the displacement variance can be obtained by expanding 
N(s)/s 2 into partial fractions and then inverting. In the case i) /3 ^ 4 , we obtain 



(X 2 (t))=2D {t-2 l(3aet 



+ a 



e 



7T 

(3 32) 

X\ [1 - E 1/2 (X-^t/o~)] -X 3 _[l- Eip jX+y/tfe" " 

(X+-X-) 




Thus, the displacement variance is proved to maintain, for sufficiently long times, 
the linear behaviour which is typical of normal diffusion (with the same diffusion coef- 
ficient as in the classical case). However, the Basset history force, which is responsible 
of the algebraic decay of the velocity correlation function, induces a retarding effect in 
the establishing of the linear behaviour of the displacement variance. As we shall see 
hereafter, this retarding effect is more evident when the Brownian particle is lighter, 
such as to give rise to regimes of effective fast anomalous diffusion characterized by 
the law 

(X 2 (t)) ~2D a t a , D a = aD(a p ) l - a - 1 0< a < 1 , Ka<2. (3.33) 



3.4 Numerical Results and Discussion 

In order to get a physical insight of the effect of the Basset history force (cou- 
pled with the added mass) on the classical Brownian motion, we exhibit the results 
obtained recently by Mainardi and Tampieri [106] concerning plots of the velocity 
autocorrelation (3.30) and the displacement variance (3.32). As an example we con- 
sider relatively light Brownian particles, by assuming x = 0.1 and x = 0.5 . We take 
non-dimensional quantities, by scaling the time with the decay constant a p of the 
classical Brownian motion and the displacement with the diffusive scale (Dap) 1 / 2 . 
Please note that here we have preferred to scale the time with o~ p more than with a e , 
since in the classical approach the added mass effect is neglected! With these scales 
the asymptotic equation for the displacement variance reads (X 2 (t) ) ~ 2t . 

In Figs 3-1 and 3-2 we plot versus the normalized time the velocity autocorrelation 
normalized with its initial value ( V 2 (0) ) and the displacement variance normalized 
with its asymptotic value 2 1 . We compare any function, provided by our full hy- 
drodynamic approach (added mass and Basset force), in continuous line, with the 
corresponding one, provided by the classical analysis, in dashed line, and by the only 
effect of the added mass, in dashed-dotted line. For large times we also exhibit the 
asymptotic estimations (3.26) and (3.31), in dotted line, in order to recognize their 
range of validity. 
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Fig. 3-1 

Velocity autocorrelation versus time for x = 0.1 (left) and for x = 0.5 (right): 
full hydrodynamic - - ; added mass classical . 
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The retarding effect of the Basset force is more evident when the Brownian particle 
is lighter, such as to appear a manifestation of fast anomalous diffusion. In fact, if 
we consider a time interval (say two decades) starting when the classical analysis 
foresees the establishment of the asymptotic linear behaviour for the displacement 
variance, a law of anomalous diffusion (X 2 (t)) ~2at a , a > 0, a / 1, can well 
approximate the exact formula (3.32), provided by the full hydrodynamic model. By 
evaluating the parameters of the anomalous diffusion, a and a , with a best fit based 
on the least squared method, we find < a < 1 and 1 < a < 2 . We recognize that 
the effective anomalous diffusion turns out to be fast; in particular, it is faster as x 
is smaller, with parameters a — > + and a — > 2~ as x — )• + . Of course, the normal 
diffusion is recovered as x — )• oo , since a — >■ l - and a — > 1 + . In Fig. 3-4 we show 
the function (X 2 )/2 versus time (in the 2-decade range 10 1 -j- 10 3 ) corresponding 
either to our analysis and to the classical analysis. While the classical curve, in 
dashed line above, is practically coincident with the linear one (regime of normal 
diffusion) our curve, in continuous line below, is fitted with a power-law curve, in 
dashed-dotted line, with an exponent a > 1 (regime of fast anomalous diffusion). 

10 3 r , , , , ;, , 10 3 r ■ ■ ■ ■ VI 




10 1 10 2 10 3 10 1 10 2 10 3 

Fig. 3-3 

The displacement variance at large times for x = 0.1 , 0.5 : 

full hydrodynamic ; best-fit — classical . 

From the above analysis we conclude that if an observer investigates the time 
evolution of a cloud of sufficiently light Brownian particles, he recognises that the 
normal diffusion is preceded by a regime of fast anomalous diffusion, which lasts for 
long time. If the observation interval is not sufficiently long, he may be induced to 
trust in the occurring of fast anomalous diffusion. 
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4. THE FRACTIONAL DIFFUSION- WAVE EQUATION 
4.1 Introduction 

By fractional diffusion-wave equation we mean the linear integro partial differen- 
tial equation obtained from the classical diffusion or wave equation by replacing the 
first- or second-order time derivative by a fractional derivative (in the Caputo sense) 
of order (3 with < (3 < 2. In our notation it reads 

B^u B^u 

— = V — , u = u(x,t), 0</3<2, X?>0, (4.1) 

where V denotes a positive constant with the dimensions L? , x and t are the 
space-time variables, and u = u(x, t) is the field variable, which is assumed to be a 
causal function of time, i.e. vanishing for t < . From the Chapter of Gorenflo and 
Mainardi [8], see in §1.3 Eq. (1.17), we remind the definition of the Caputo fractional 
derivative of order (3 > for a (sufficiently well-behaved) causal function f{t) , 



i r* 

T(m - P) J 



D?f(t) := J m ^D m f(t) = 

( dt m 

Introducing the causal power function 



(t- T y+ i - m 

d m 



dr , m — 1 < f3 < m , 



f(t), p = m. 



ffA)' A>0 ' 

where the suffix + is just denoting that the function is vanishing for t < , and 
recalling the Laplace transform pair $>\(t) -j- s~ x , we easily recognize that 

m— 1 

DZ f(t) := <S> m - P {t) * /<">(*) -h J f(s) - s P ~ 1 - k f (k) (0+) , m - 1< < m . 

We note $>\(t) * $ M (t) = Qx+^t) . In Eq. (4.1) we thus need to distinguish two cases 
i) < j3 < 1 , and ii) 1 < f3 < 2 , for which the equation assumes the explicit 
forms as follows : 

The equations (4.2) and (4.3) can be properly referred to as the time-fractional dif- 
fusion and the time-fractional wave equation, respectively. 
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A fractional diffusion equation akin to (4.2) has been explicitly introduced in 
physics by Nigmatullin [110] to describe diffusion in special types of porous media, 
which exhibit a fractal geometry. The author [111] has shown that the fractional 
wave equation governs the propagation of mechanical diffusive waves in viscoelastic 
media which exhibit a simple power-law creep. This problem of dynamic viscoelas- 
ticity, formerly treated by Pipkin [5] but unaware of the interpretation by fractional 
calculus, thus provides an interesting example of the relevance of (4.3) in physics. 
Of course, anytime some hereditary mechanisms of power-law type are present in 
diffusion or wave phenomena, the appearance of time fractional derivatives in the 
evolution equations is expected. 

In a series of papers [112-116] the author has pursued his analysis on the fractional 
diffusion- wave equation (4.1), based on Laplace transforms and special functions of 
Wright type. Mathematical aspects of integro differential equations akin to (4.2-3) 
and based on the use of integral transforms and special functions have been treated 
in some relevant papers by Wyss [117], Schneider and Wyss [118], Schneider [119] 
(Mellin transforms and Fox H functions) and by Fujita [120] (Fourier transforms 
and Mittag-Leffler functions). More formal approaches based on semigroup theory in 
Banach spaces have been given by Kochubei [121-122] and El-Sayed [123]. Recently 
the integro-differential equation treated by Fujita has been considered by Engler [124] 
in a very interesting paper in view of the connection between similarity solutions and 
stable probability distributions. 

Hereafter we present a review on the fractional evolution equation (4.1), essentially 
based on our works [111-116]. In §4.2 we analyse the two basic boundary- value 
problems, referred to as the Cauchy problem and the Signalling problem, by the 
technique of the Laplace transform and we derive the transform expressions of the 
respective fundamental solutions (the Green functions). 

In §4.3 we carry out the inversion of the relevant Laplace transforms and we 
outline a reciprocity relation between the Green functions themselves in the space- 
time domain. In view of this relation the Green functions can be expressed in terms 
of two interrelated auxiliary functions in the similarity variable r = \x\/ {\fDt^' 2 ) . 
These auxiliary functions can be analytically continued in the whole complex plane 
as entire functions of Wright type. 

In §4.4 we show the evolution of the fundamental solutions of both the Cauchy 
and Signalling problems for some (rational) values of the order of time derivation. 
To gain more insight into the phenomenon of fractional diffusion we also exhibit the 
evolution of an initial box function in the Cauchy problem. This allows us to better 
recognize the processes of slow diffusion (0 < f3 < 1) and the intermediate processes 
between diffusion and wave propagation (1 < j3 < 2). 
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Finally, in the Appendix, we provide the reader with a review of the main mathe- 
matical properties of our auxiliary functions in the framework of the Wright functions. 
The interesting connection with the stable probability distributions is not treated but 
deferred to recent works of ours, there quoted. 

4.2 Analysis of the Cauchy and Signalling Problems with the Laplace Transform 

As well known, the two basic boundary-value problems for the evolution equations 
of diffusion and wave type are the Cauchy and Signalling problems. In the Cauchy 
problem, which concerns the space-time domain — oo < x < +00 , t > , the data are 
assigned at t = + on the whole space axis (initial data). In the Signalling problem, 
which concerns the space-time domain x > , t > , the data are assigned both 
at t = + on the semi-infinite space axis x > (initial data) and at x = + on 
the semi-infinite time axis t > (boundary data); here, as mostly usual, the initial 
data are assumed to be vanishing. Extending the classical analysis to our fractional 
equation (4.1), and denoting by f(x) and h(t) two given, sufficiently well-behaved 
functions, the basic problems are thus formulated as following: 

a) Cauchy problem 

u(x,0 + ) = f(x) , -oo<£<+oo; w(=Foo,t) = , t>0; (4.4a) 



b) Signalling problem 

u(x,0 + ) = 0, x>0; 



u 



(0+,t) = h(t), u(+oo,t) = 0, £>0. (4.46) 



If 1 < f3 < 2 , we would add in (4.4a) and (4.4b) the initial value of the first-order 
time derivative of the field variable, i.e. -^:u(x,0 + ) = g(x) , since in this case Eq. 
(4.1) turns out to be of the second order in time, see the integro-differential equation 
(4.3), and, consequently, two linearly independent solutions are to be determined. 
We limit ourselves to choose g(x) = . We easily recognize that the above Cauchy 
problem for (4.1) can be expressed through the integral equations of fractional order 



f(x) + 



V 



u(x,t) = < 



f(x)+tg(x) + 



d 2 u 
dx 2 
V 



(x,t) 



< < 1, 



d 2 u 
dx 2 



(x,t) 



(4.5) 



{t-rf^dT, l</3<2. 



We thus note that for 1 < (3 < 2 the choice g(x) = ensures the continuous de- 
pendence of the solution on the parameter (3 also in the transition from (3 = 1~ to 

In view of our analysis we find it convenient to put 



v = ^ , < v < 1. 



(4.6) 
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For the Cauchy and Signalling problems we introduce the so-called Green func- 
tions Q c (x, t; v) and Q s (x, t; v), which represent the respective fundamental solutions, 
obtained when g{x) = 8{x) and h(t) = S(t) . As a consequence, the solutions of the 
two basic problems are obtained by a space or time convolution according to 

/+oo 
Gc{x-tt;v)f{g)<%, (4.7a) 
-oo 

f t+ 

u(x,t;u) = / Q s {x,t — t;u) h(r) dr . (4.76) 
Jo- 
lt should be noted that Q c (x,t; v) = Q c (\x\, t; v) since the Green function turns out 



to be an even function of x. 

For the standard diffusion equation {y = 1/2) it is well known that 

&(M;l/2) :=£ c d (M) = _^t-V2 e -x 2 /(4in) ^ (4ga) 

2v tiT) 

1/2) :=Qt(x,t) = _J^ t -3/i e -x 2 /(4Vt) 
2v tvD 

For the standard wave equation (v = 1) it is well known that, putting c = \fV> , 
g c {x,t;l) :=g™(x,t) = ^ [6{x - ct) + 5{x + ct)} , (4.9a) 

g s (x,fl) :=g™(x,t)=5(t-x/c). (4.96) 



In the general case < v < 1 the two Green functions will be determined by using 
the technique of the Laplace transform. This technique allows us to obtain the trans- 
formed functions G c (x, s; u), g s (x, s; u), by solving ordinary differential equations of 
the 2-nd order in x and then, by inversion, g c (x, t; v) and g s (x, t; v). 

For the Cauchy Problem (4.4a) the application of the Laplace transform to (4.1) 
with u(x,t) = g c (x,t;u), and g c (x,0 + ; v) = f(x) = 5(x) , [and -^g c (x, + ; v) = 
\{ 1/2 < v <1] leads to the non-homogeneous differential equation satisfied by the 
image of the Green function, G c (x, s; v) , 

V 4% - s^gl = -S(x) s 2v - 1 , -oo < x < +oo . (4.10) 
ax 2 

Because of the singular term 5(x) we have to consider the above equation separately 
in the two intervals x < and x > 0, imposing the boundary conditions at x = =Foo , 
(? c (=Foo,t; v) = , and the necessary matching conditions at x = ± . 
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We obtain 

(T c (x, s; v) = e~{\ x \/^V -oo<x<+oo. (4.11) 

2VVS 1 -" 

In fact, from (4.10) Q c is expected in the form 

_ ( Cl (s)e-W^ sV +C2 (s)e+^/^ s \ if x>0; 

G c (x,s;v)= I _ (4.12) 

[c 3 (s)e-W^ sV + C4 (s)e+^/^ s \ if x<0. 

Clearly, we must set C2(s) = C3(s) = , in order to ensure that the solution vanishes 
as |x| — > oo . We recognize from (4.10) that in x = the function Q c (x,s;v) is 
continuous but not its first derivative: we write 

&(0+ s; z/) - K(0-,s; i/) = ci(s) - c 4 (s) = , (4.13) 

and, by integrating (4.10) with respect to a; from a; = 0~ to x = + , 

^ &(0+ S ; i/) - A &( - s; !/) = _[ Cl ( a ) + C4 ( S )] = • (4.14) 

Therefore, using (4.13-14) we obtain ci(s) = c 4 (s) = l/(2 v / ^s 1 " i ') , and conse- 
quently the expression (4.11). 

For the Signalling Problem (4.4b) the application of the Laplace transform to (4.1) 
leads to the homogeneous differential equation 



with u(x,t) = G s {x,t;v) , g s (x,0+;v) = 0, [and 4£ s (a;, 0+; v) = if 1/2 < v < 1 



rl 2 G ~ 

P^-s 2 ^ s = 0, x>0. (4.15) 

Imposing the boundary conditions at x = 0, £ s (0 + ,t;z/) = h(i) = 8(t) , and at 
x = +oo , Q s (+oo, £; z/) = , we obtain 

&(x,s;i/) = e~( x / v ^) s '\ x>0. (4.16) 

In fact, from (4.15) Q s is expected in the form 

(T a (x, s; v) = ci(s) e"^^ s " + c 2 (s) e+^Z^ s \ x > . (4.17) 

Clearly, we must set C2(s) = to ensure that the solution vanishes as x — > +oo , and 
consequently we obtain c\(s) = Q s (0, s; v) = S(s) = 1 . 
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4-3 The Reciprocity Relation and the Auxiliary Functions 
From (4.11) and (4.16) we recognize 

■%-(Ts = -2vx'g c , £>0, (4.18) 
as 

which implies for the original Green functions the following reciprocity relation 

2uxQ c (x,t;u) = tQ a (x,t;u), x>0, t>0. (4.19) 



The above relation can be easily verified in the case of standard diffusion {y = 1/2) 
where the explicit expressions (4.8) of the Green functions leads to the identity (for 
x > 0, t > 0) 

xg*{x,t) = tgl{x,t) = 7 ^ r ^= e- x2 /( AVt ) =F d (r) = - M d (r) , (4.20) 
where 

r = x/(\ r Dt 1/2 ) > 0, (4.21) 
is the well-known similarity variable and 

M d (r) = -^e- r2 / 4 . (4.22) 
v 7r 

We can refer to F d {r) and M d (r) as to the auxiliary functions for the diffusion 
equation because each of them provides the fundamental solutions through (4.20). 
We note that M d (r) satisfies the normalization condition J °°M d (r) dr = 1 . 

Now we are going to show how, in the general case < v < 1 , the inversion of the 
Laplace transform in (4.11) or (4.16) leads us to generalize the auxiliary functions 
F d (r) and M d (r) by introducing the proper similarity variable for x > , t > , 

r = x/(Vvt u )>0. (4.23) 



The new auxiliary functions, that we denote by F(r; v) and M(r; v) , turn out to 
be expressed in terms of Bromwich complex integrals as shown hereafter. 

Applying in the reciprocity relation (4.19) the complex inversion formulas for the 
transformed Green functions (4.11) and (4.16), we obtain 

2uxg c (x,t;u) = ^=^- [ e st-(x/Vv)s»^ X>Q f>Q (424) 
V ' 2VV 2ttz J Br s 1 -" v ' 
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and 



tQ s {x,t-u) = t— I e st ~ ( x /^) s " ds , x>0, *>0, (4.25) 
2m J Br 

where Br denotes the Bromwich path. 

In order to express the Bromwich integrals in terms of the single similarity variable 
r defined by (4.23) let us change the integration variable in (4.24-25) setting a = st . 
From (4.24) we obtain 

2vxG c (x,t;v) = vrM(r;v) , (4.26) 



'{r;v):=±, / r>0, <i/<l, (4.27) 

^1 J Br G 



with 

' Br 

and, from (4.25), 

tg 8 (x,t;u) = F(r;v), (4.28) 

with 

F(r;v):=— [ e a ~ ra " do , r>0, < v < 1 . (4.29) 
2ttz J Br 

Therefore we conclude that, for x>0,t>0,r>0, 

2vxQ c {x,t\v) =tg s (x,t;u) = F(r; v) = vrM(r;v). (4.30) 

The above definitions of F(r; v) and M(r; z/) by the Bromwich representation can 
be analytically continued from r > to any z G C, by adopting suitable integral and 
series representations valid in all of C. 

For this purpose, let us deform the Bromwich path Br into the Hankel path Ha 
[a contour consisting of pieces of the two rays argcr = ±<p extending to infinity, 
and of the circular arc a = ee* 9 , \6\ < <f> , with G (n/2,n)] which is chosen to 
be equivalent to the original path (at least) for z real and positive. The Hankel 
integral representation allows us to obtain the series representation for each auxiliary 
function. In fact, after expanding in series of positive powers of z the exponential 
function, exp(— za") , exchanging the order between the series and the integral and 
using the Hankel representation of the reciprocal of the Gamma function, 



7n=^- [ (gC, 
(0 2m J Ha 



r(C) 

we finally obtain the required power series representation. Since the radius of con- 
vergence of the power series can be proven to be infinite for < v < 1, our auxiliary 
functions turn out to be entire in z and therefore the exchange between the series 
and the integral is legitimate. 
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The integral and series representations of F(z; u) and M(z; u), valid on all of C , 
with < v < 1 turn out to be 



' 1 

2tt7 



do 



F(z;u) = { 



oo 



E 



Ha 

' nl T(—vn) 



zGC, < v < 1 , 



(4.31) 



and 



' 1 

2ni 



a-za v d ° 



l-u 



M(z;v) = < 



E 

n=0 



Ha ° 

n\T[-vn+ (1 - u)] 



zeC, 0<z/<l 



(4.32) 



In Appendix we show that the two auxiliary functions turn out to be particular 
examples of a special entire function known as Wright function. We refer the reader 
to the Appendix for the main properties of the auxiliary functions including the 
related Laplace transform pairs. 

44 Plots 

In the following figures we exhibit the plots of the Green functions for both the 
Cauchy and Signalling problems, in the cases f3 = 1/2, 1 , 3/2, (j3 = 2u) takining 
V= 1. 

The plots of G c {x,i) versus x at fixed times (t = 1,2,3) for (3 = 1/2, 1,3/2 are 
reported in Figs 4-1, 4-2, 4-3, respectively. The plots of Q s {x,i) versus t at fixed 
positions (a; = 0.9,1,1.1) for (3 = 1/2,1,3/2 are reported in Figs 4-4, 4-5, 4-6, 
respectively. 
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G c (x,t) 


t=1 




P=3/2 






\ t=2 










x 



12 3 4 



Fig. 4-3. 

The Green function Q c (x, t) for f3 = 3/2 versus x at fixed times. (Cauchy problem) 



G s (x,t) 




\x=0.9 


P=1/2 




t 



0.02 0.04 0.06 0.08 0.1 



Fig. 4-4 

The Green function Q s (x,t) for f3 = 1/2 versuts t at fixed positions. (Signalling problem) 
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Fig. 4-5 

The Green function Q s (x, t) for j3 = 1 versuts t at fixed positions. (Signalling problem) 
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0.5 



G s (x,t) 




x=0.9 






(3=3/2 


/ /x/f\. 






P 



0.5 1 1.5 2 



Fig. 4-6 

The Green function Q s (x,t) for f3 = 3/2 versuts t at fixed positions. (Signalling problem) 

In order to gain more insight into the phenomena governed by the fractional 
diffusion wave equation (4.1), we consider a simple Cauchy problem where the initial 
data are provided by a box-type function. Precisely, taling V = 1 , we consider 

r+l 

u(x,0 + ) = 6(1 - u(x,t) = J g c (x-£,t;v)d£. (4.33) 

In the following we exhibit plots of the solution versus x (0 < x < 3), at fixed 
t (t = , 0.5, 1), for some fractional values of /3 = 2v . In Fig. 4-7 we compare the 
cases concerning the fractional diffusion equation (0 < /3 < 1), whereas in Fig. 4-8we 
compare those concerning the fractional wave equation (1 < (3 < 2). 
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u(x,t) 


t=1 







12 3 




(3=2/3 u ^ 



0.5 



12 3 



12 3 




12 3 




3 12 3 
Fig. 4-7 

Evolution of the initial box-signal (dashed line) at t = 0.5 (left) and t = 1 (right), 
versus x , for various values of/3: 1/2,2/3,1. 

In Fig. 4-7 we thus obtain some comparison between the fractional diffusion and 
the standard diffusion. We easily recognize for < (3 < 1 a diffusive behaviour, which 
is slower with respect to the case /3 = 1 of standard diffusion: this is consistent with 
a slow diffusion process. 
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3 

Fig. 4-8 

Evolution of the initial box-signal (dashed line) at t = 0.5 (left) and t = 1 (right), 
versus x , for various values of/3: 4/3,3/2,2. 

In Fig. 4-8 we recognize the intermediate process between standard diffusion 
(where discontinuities are smoothed out) and wave propagation (where discontinuities 
can propagate with finite speed). This can be seen from the appearance of a hump, 
which tends to be narrower as (3 — > 2~ up to reproduce the discontinuities of the 
signal for (3 = 2. For 1 < (3 < 2 the hump travels with finite velocity (as in a wave 
process) but the signal diffuses instantaneously (as in a diffusion process). 
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APPENDIX: THE WRIGHT FUNCTION 

In this Appendix we shall consider the general class of the Wright functions with 
special regard to the special functions that we have introduced for the sake of con- 
venience in the treatment of the fractional diffusion wave equation, the so-called 
auxiliary functions. It is our purpose to provide a review of the main properties 
of these functions including their series and integral representations and the related 
Laplace transform pairs. We also mention their connection with the Mittag-Leffler 
functions, for which we refer the reader to Gorenflo and Mainardi [8]. 

A.l The representations for the Wright function W\ jfl (z) 

The Wright function, that we denote as W\^(z) , where A > — 1 and fx > , is so 
named from the British mathematician E. Maitland Wright who in 1933 introduced 
it in the asymptotic theory of partitions. A list of formulas concerning this function 
can be found in the handbook of the Bateman Project [125]. We note that originally 
Wright considered A > [126] and only later, in 1940, he extended to — 1 < A < 
[127]. Relevant investigations on this functions have been carried out by Stankovic 
[128-129], by other authors quoted by Kiryakova [130], and more recently by Gorenflo, 
Luchko and Mainardi [142-143] and by Wong and Zhao [144-145]. 

The Wright function is defined by the series representation, valid in all of C , 



zgC, A > -1 , fx>0, 



(A.l) 



so that it turns out to be an entire function. This property remains valid even if \x 
is an arbitrary complex number. The integral representation reads 

W x Jz) = — [ e ^ + za- x ^ zeC A> _! ^ >0 ( A2 ) 
,M J 2ni J Ha a* J 

where Ha denotes the Hankel path. The formal equivalence between the two repre- 
sentations is easily proved using the Hankel formula for the Gamma function 

1 

'Ha 



no 



e u u c du. 



and performing a term-by-term integration. In fact 

-a da 1 



oo 

= E 



+ zo~~ 



0-M 



2ni 



Ha 



T - 

,n=0 



a 



da 
a^ 



n=0 



Z' b 

nl 



1 

2ni 



e a a- Xn -> 1 da 



Ha 
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It is possible to prove that the Wright function is entire of order 1/(1 + A) , hence 
of exponential type if A > . The case A = is trivial since W ^(z) = e z /T(ji) . 

Wright showed in particular that in the case A = — v G (—1,0) there is the 
following asymptotic expansion, valid in a suitable sector about the negative real 
axis 

W-^{z)=Y 1 ' 2 -»e- Y ( A m Y- m + 0(\Y\- M )\ as z -+ -oo , (A3) 

\m=0 / 

with Y = Y{z) = (1 — u) {—v v z) 1 ^ 1- ^ , where the A m are certain real numbers. 

A. 2 The Wright functions as generalization of the Bessel functions 

For A = 1 and n = v + 1 the Wright function turns out to be related to the well 
known Bessel functions J u and I u by the following identity 

(z/2y Wh „ +1 ^zy4) = { J ; {z) . (A4) 

[ w) 

In view of this property some authors refer to the Wright function as the Wright 
generalized Bessel function (misnamed also as the Bessel-Maitland function) and 
introduce the notation 

n=0 V ' 

As a matter of fact, the Wright function appears as the natural generalization of the 
entire function known as Bessel - Clifford function, see e.g. [130], and referred by 
Tricomi [131-132] as the uniform Bessel function *, 

T(z; v) := z~^ J v {2jz) = £ „, r ^ { * = W^ +1 {-z) . (A6) 

n=0 

Some of the properties which the Wright functions share with the most popular 
Bessel functions were enumerated by Wright himself. Hereafter, we quote some rel- 
evant relations from the Bateman Project [125], which can easily be derived from 
(A.l): 

XzW x ,x+M = Wx,^-i(z) + (1 - n) W X ,M > (A-7) 
^W x , fl (z)=W x , x+fl (z). (A8) 



* The great Italian mathematician denoted this function by E v {z) . Here we write T{z\ v) to 
remain in accordance with the standard notation used for the Mittag-LefHer function [8] . 
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A. 3 The Auxiliary Functions F(z; v) and M(z; v) 

In our treatment of the time fractional diffusion wave equation we have found it 
convenient to introduce two auxiliary functions F(z; v) and M(z; v) , where z is a 
complex variable and v a real parameter < v < 1 . Both functions turn out to be 
analytic in the whole complex plane, i.e. they are entire functions. Their respective 
integral representations read, see (4.31-32), 

F( z -v):=^- [ e° za " da , zGC, < v < 1 , (A9) 



2ni 



Ha 



M(z-u):=— [ e a ~ zaU , zeC, < v < 1 . (A10) 

From a comparison of (A. 9-10) with (A. 2) we easily recognize that these functions 
are special cases of the Wright function according to 

F{z;v) = W- vfi {-z), (All) 

and 

M(z;v) = W_„, 1 -„(-z). (A12) 
From (A. 7) and (A. 11-12) we find the relation 

F(z\v) = vzM(z\v). (A13) 

This relation can be obtained directly from (A. 9-10) with an integration by parts, 
i.e. 

e° = I e" f-J-Ae--^ da = -1 f e° - ^ 'da . 



Ha O- 1 v J H ~ \ VZ da~ J " VZ J Ha 

The series representations for our auxiliary functions turn out to be respectively, 
see (4.31-32), 

00 (—z) n 1 °° (— ?) n 

F & ■= E ^W^y = E v{un + 1} sin(7n/n) ' (A14) 

n=l ^ ' n=l 

and 

M(z; u) := V —. [ ~ Z \ ^ = ~Y r, T{yn) sin(7rz/n) . (A15) 

n=0 L V n n=l v ; 

The series at the R.H.S. have been obtained by using the well-known reflection for- 
mula for the Gamma function T(Q F(l — Q = tt/ sin n( . Furthermore we note that 
F(0; v) = , M(0, v) = 1/T(1 — v) and that the relation (A. 13) can be derived also 
from (A.14-15). 
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Explicit expressions of F(z; v) and M(z; v) in terms of known functions are ex- 
pected for some particular values of v. Mainardi & Tomirotti [114] have shown that 
for v = 1/q, where q > 2 is a positive integer, the auxiliary functions can be ex- 
pressed as a sum of (q — 1) simpler entire functions. In the particular cases q = 2 
and q = 3 we find 

M fel /2) = £(-1)- -) = ^exp(-V4) , (A!6) 

v m=0 m 

and 

1 °° /l\ 2 3m 1 °° /2\ ^Sm+l 

M(z; 1/3) = WN) £ U J m Jsmy. ~ Wm £ U J m J^TTy. (A17) 

^S^Ai^ 1 / 3 ) , 

where Ai denotes the Airy function. 

Furthermore it can be proved [114] that M(z; 1/q) satisfies the differential equa- 
tion of order q — 1 

di- 1 (-l) q 

M(z; 1/q) + z M(z; 1/q) = , (A18) 

subjected to the q — 1 initial conditions at z = 0, derived from (A. 15), 

M W (0;l/g) = fclL r[(/i + l)/g] sm[n(h + l)/q], h = 0, 1, . . . q - 2 . (A19) 

We note that, for q > 4 , Eq. (A. 18) is akin to the hyper-Airy differential equation 
of order q — 1 , see e.g. [133]. Consequently, in view of the above considerations, the 
auxiliary function M(z; v) can be referred to as the generalized hyper- Airy function. 

Let us now consider the problem of the asymptotic evaluation of the function 
M(z; v) as \z\ — > oo in the complex plane. Referring to a preliminary report of ours 
[134] for the detailed asymptotic analysis in the whole complex plane, which includes 
the phenomenon of Stokes lines, here we limit ourselves to provide the asymptotic 
representation as z = r is real and positive by using the ordinary saddle-point method. 
Choosing as a variable r/v rather than r the computation is easier and yields, see 
[H4], 

M{r/v- v) ~ a(u) A u ~ 1 / 2 )/ ^ ~ v ) exp \-b(u) rV(l ~ r -+ +oo , (A20) 

where 

a(u) = -7=== > , b{u) = — ->0. (A.21) 
a/2/T (1 — v) v 

The above evaluation is consistent with the first term in Wright's asymptotic expan- 
sion (A. 3) after having used (A. 12). 
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The exponential decay for r — > +00 ensures that M(r; v) is absolutely integrable 
in M + . We can easily prove the normalization property in M + 

r+oo 

/ M(r; v)dr = l, (A.22) 
Jo 

and more generally we can compute all the moments in 1R + 



L 



+ °° Yin + 1) 

. r n M{r-u)dr= / \ , n e N . (A.23) 

The results (A. 22-23) can be formally derived by using the Laplace transform of 
M(r; z/) , as shown later. Analogously we can compute all the moments in M + for 
F(r; v) . 

A. 4 The Laplace transform pairs related to the Wright function 

Let us consider the Laplace transform of the Wright function using the following 
notation 

/■OO 

WxA±r) + C [W x ^{±r)\ := / e ~ s r W x ,»{±r) dr , 

where r denotes a non negative real variable, i.e. < r < +00 , and s is the Laplace 
complex parameter. 

When A > the series representation of Wright function can be transformed 
term-by-term. In fact, for a known theorem of the theory of the Laplace transforms, 
see e.g. Doetsch [109], the Laplace transform of an entire function of exponential type 
can be obtained by transforming term-by-term the Taylor expansion of the original 
function around the origin. In this case the resulting Laplace transform turns out 
to be analytic and vanishing at infinity. As a consequence we obtain the Laplace 
transform pair 

W x ^(±r) - - s E x ^ (±^J , A>0, |s|>p>0, (A24) 

where E x ^ denotes the generalized Mittag-Leffler function in two parameters, and p 
is an arbitrary positive number. The proof is straightforward noting that 

^ n\ Y{\n + p) • ~s ^ Q F(\n + p) ' 
and recalling the series representation of the generalized Mittag-Leffler function, 

E^iz)-^^— — , a>0,i/eC, zeC. 

r (an + v) 

n=0 V ' 
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For A — > + (A. 24) we obtain the Laplace transform pair 



±r 



1 



1 



1 



E, 



0,fJ, 



±- 



1 



r(/x) r(//) s =f i s 

where, to remain in agreement with (A. 24), we have formally put 

1 _ , . 1 1 



s\ > 1 , 



(A25) 



oo 



n=0 vry 



r(/*) 



^o(^) 



r(/i) i-z 



We recognize that in this special case the Laplace transform exhibits a simple pole 
at s = ±1 while for A > it exhibits an essential singularity at s = . 

For — 1 < A < the Wright function turns out to be an entire function of order 
greater than 1, so that care is required in establishing the existence of its Laplace 
transform, which necessarily must tend to zero as s — > oo in its half-plane of conver- 
gence. For the sake of convenience we limit ourselves to derive the Laplace transform 
for the special case of M(r; v) ; the exponential decay as r — > oo of the original 
function provided by (A. 20) ensures the existence of the image function. From the 
integral representation (A. 10) we obtain 



M{r;v) -=- 



1 

2tH 
1 

2ni 



-s r 



Ha 



C 
Ha 

oo 



a 



ro 



v da 



l-v 



a 

-r{s + or v ) 



dr 



dr 



1 r e a a' y - 1 
2^1 J Ha a- + s 



da 



da . 



Then, by recalling the integral representation of the Mittag-Leffler function, 



E, 



< {z) = h L 



^a-l e C 



d( , a > , z E C , 



(<* -z 

we obtain the Laplace transform pair 

M(r;u) + E u (-s) , < v < 1. 



(A26) 



In this case, transforming term-by-term the Taylor series of M(r; v) yields a series 
of negative powers of s, which represents the asymptotic expansion of E v {—s) as 
s — > oo in a sector around the positive real axis. We note that (A. 26) contains the 
well-known Laplace transform pair, see e.g. [109], 

M(r; 1/2) := — exp (- r 2 /A) + E 1/2 (-s) := exp (s 2 ) erfc (s) , s G C . 



We also note that (A. 26) allows us to derive (A. 22-23) by accounting for the property 



+0O 



r n M(r;u)dr = lim (-1) 7 

s— >-0 



d 1 



ds 1 



E v {-s) . 
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Analogously, because of (A. 3), we can prove that in the case A = — v e (—1,0) 
we get 

W-^(-r) + E^ +u (-s) , < v < 1 . (A27) 
In the limit as A — > 0~ we formally obtain the Laplace transform pair 

'W-O-g^j^ (A28) 

Therefore, as A — > ± , we note a sort of continuity in the formal results (A. 25) and 
(A.28) since 

1 f (l/s)E (-l/s), \s\>l; 



(* + !) [ ^o(-s), |s|<l 



(A29) 



A quite relevant Laplace transform pair related to the auxiliary functions of ar- 
gument r~ v is 

-F(l/r";u) = -^M{l/r u ;v) exp (-s v ) , < v < 1 . (A30) 

We recall that a rigorous proof of (A. 30) was formerly given by Pollard [135], based 
on a formal result by Humbert [136]. The Laplace transform pair was also obtained 
by Mikusihski [137] and, albeit unaware of the previous results by Buchen & Mainardi 
[138] in a formal way. 

After noting that the pair (A. 30) can be easily deduced (with r = t) from (4.16) 
and (4.30), hereafter we like to provide two independent proofs carrying out the 
inversion of exp(— s v ) , either by the complex integral formula or by the formal series 
method. We obtain 



C- 1 [exp(-0] = — / e sr ~ s "ds = — [ e a 
1 K h 2vu J Ha 2m r J Ha 



da 



'Ha z ' nl1 JHa 

*F(l/rV) = ^M(l/rV), 



and 



n=0 n=l v ' 

= iF(l/rV) = ^M(l/rV). 

Last, but not the least, we would like to mention the relevance of our auxiliary 
functions in probability theory. In fact, as shown by Engler [124], they turn out to 
be related with the probability density functions of the so-called stable distributions. 
For this interesting topic we refer the reader to our recent works [139-141]. 
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